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ABSTRACT 


The project initiated a study of a mathematical model of a tunable 
Ti:Sappnire solid-state laser. A general matnematical model was developed 
for the purpose of identifying design parameters which will optimize the 
system, and serve as a useful predictor of the system's oenavior. 
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MATHEMATICAL MODELING OF A TI -.SAPPHIRE SOLID-STATE LASER 


By 

John J. Swetits* 

INTRODUCTION 

Recently, several new experimental techniques involving the transient 
behavior of solid state lasers have been developed to control the frequency 
and shape of a laser beam (1-4). Among these is injection seeding, a tech- 
nique which relies upon precise timing between the pumping beam, the injec- 
tion beam and the onset of lasing action. Its effective use requires accu- 
rate understanding of the temporal development of the laser dynamical vari- 
ables. It follows that there is a need to reexamine the transient behavior 
of the laser under the conditions which prevail in this situation. In order 
to study the use of injection seeding and other similar techniques, we have 
developed a model of tne dynamical process in a four-level solid-state las- 
er. The analysis was developed for the Titanium- doped Sapphire system, but 
it is of a general character, and can be applied broadly witn minor modifi- 
cation to any system. 

In solid-state lasers the optically active ions are fixed in place 

3+ 

within a transparent host. In Ti:Sapphire, the Ti ion substitutes for the 
3+ 

A1 ion in Al^O^ to a concentration of about 1%. Because of a strong coup- 
ling between the active ion and the host lattice the absorption and emission 
lines are broadened and the material can lase at over a range of wave- 
lengths. If left running freely, the laser will operate at a wavelength 
near to the peak of the emission cross-section. By introducing an injection 
signal the laser can be made to lase at a different wavelength. The 
solid-state laser material is formed into a cylindrical rod with its axis 

♦Professor, Department of Mathematical Sciences, Old Dominion University, 
Norfolk, Virginia 23508. 


coincident with the optical axis of a Fabry-Perot cavity. This rod is 
presumed to be pumped axially by another laser at the pumping wavelength 
although this mode of pumping is not critical to the model. In addition, an 
injection signal, at the desired lasing wavelength, is injected axially into 
the materi al . 

The model describes the time evolution of the occupancy of the 
electronic levels and of the lasing photon flux concentration within the 
cavity. Most of the parameters used in the model can be obtained by direct 
measurement; the effects of the few quantities which cannot be directly 
measured are studied parametrically. The main purpose of this effort is to 
understand the temporal characteristics of injection seeding on solid-state 
lasers and so the model avoids some complications of laser operation wnicn 
are not essential to this objective. For example, it treats a single lasing 
wavelength with a narrow bandwidth. It also does not directly treat the 
spatial distribution of excitation in the rod, but considers the spatial 
average over the rod as a whole. 

In the next section, the basic system of equations for the electron and 
photon populations are introduced and put into the form finally used. 
Following that, the quantitative properties of the solutions to this system 
is derived. The system is proven to be asymptotically stable and general 
bounds on tne population inversion are obtained. After that a nunerical 
analysis of the system is described and solutions obtained. 

THE RATE EQUATIONS 

The emission and absorption spectra of Ti 3+ :Al 2 0j (5-6) are the basis 
for an idealized four-level model of the lasing action of Titaniim in 

Sapphire shown in Fig. 1. 
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Figure 1. Energy level diagram for the model four level laser. 

Only those non-radiati ve transition to or from the 
lasing levels are shown in this figure. 
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The rate equations giving the time dependence of the population 
densities, n^t), for each of the electronic levels for this system are 
given below. 
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In these equations is the lifetime for the i electronic level, 
while is the rate for the transition from level i to level j. 

1 J 

These rates are related by reciprocal sum rule 


W 0 3 ( t ) the Pimping rate giving the number of photons per second pumped to 

t h 

level 3 from the ground state. The degeneracy of the i tn electron level is 
g l - * The population densities, n. , for each of the laser levels are 
constrained by the relation 


n T0TAL 


VVW 


(4) 
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wnere n jQ-p^|_ is the dopant concentration. The photon concentration 
(number of photons per cm ) within tne lasing medium is represented Dy 
the emission cross-section by a and the speed of light in the crystal by v. 
The rate of change of the photon concentration is given oy the following 
expression: 




n^) a<f>v - — + S f 


(5) 


Here t is the cavity lifetime and it represents the loss of photons from 
c 

all mechanisms in the cavity; S is the contribution to the photon flux in 

P 

the laser beam coming from the spontaneous emission of excited atoms. It is 
expressed as a fraction Spo of the total fluorescence decay rate, . Spo 
depends both on the geometry of the laser cavity and on the fluorescence 
spectrum of tne lasing ion. 

We can simplify these equations somewhat without losing generality. 
First, we assume that the non-radi ative decay from the pumped level to the 
upper lasing level is fast, t^>> all other time scales. Then, integrating 
the rate equation for n with this assumption, we obtain 

n 3 U) - W 0 J (t)n 0 (t)T 3 


so that n 3 can be replaced by this expression in the subsequent equations. 
This condition is equivalent to assuming that dn^/dt=0. Furthermore, we 
assume that the dominant decay path from level 3 is the non-radi ative decay 
to level 2, that is t ^ <<t.^ and T 3 g» so that 1° this case tne 

pumped photons are transferred to the upper laser level at a rate 
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(n 3 / T 3 2 ) " ^O3 n 0*- T 3 / ' T 32'* W p n 0’ 


and n /t = 0. Furthermore, if n <<n. for i = 0, 1, or 2, so that not many 
electrons stay in their pumped level 

n = n +n +n +n * n +n + n • 

TOTAL 0 1 2 3 0 2 2 

Now the decay from level 2 has both radiative and non-radiative 
contributions wnich are related by 

t 2 1 = T fl 1+t NR 1 * ^ 

With these assumptions the rate equations for the electron concentrations 
become: 


dn 2 g 


dt 


n Q - (n 2 - A n x ) a<f>v - J. 
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t 2 


(7) 


dn. n_ g. 

= — + (n 2 n x ) a<t>v - _ 


dt t 


fl 
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( 8 ) 


n T0TAL = n Q + n i + n 2 


(9) 


Equations (5), (7), (8) and (9) form the basis of this model. We reorganize 
tnem sligntly. First, we introduce new variables replacing the laser level 
populations, n^ and n 2> with the population inversion n = n^g^g^n^, and 
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n l» th e lower laser level population. Next we scale tnese new variables and 
the photon concentration to tne total cocentration of ions, n jQj i t\L 

x * n/n T0TAL} 

y = n l/ n T0TAL, 3nd 

2 = * /n T0TAL - 

Finally, we scale the time variable to the time interval, x^, taken 
here to be one nanosecond. With these changes the laser dynamic equations 
assume the form 


• U 

v = Cv + a(-y) (xz) + ( p (l-x-Yy)+ 
1 0 
1 I(t) 


In this expression, v is a vector with components equal to the three 
normalized concentrations (x,y,z). The matrix C representing tne linear 
part of the differential equation has components given below. 


C . - (-L + L±) c 12 = (l-Y) c 13 - 0 

t 2 T f* t 2 t U t 1 


•21 


'31 



1 

T-l 1 

= 


C 22 7“ “ 7“ 


T ft 

T ft T 1 

— 

Spo 

c 32 = ^ y ‘ 1 ^ C 31 


T f£ 



C 23 = 0 


'32 
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The second term has the non-linearity of the system as the product of the 
population inversion and photon concentration. The strength of this non- 
linearity is given by the parameter. 

° ’ 0 v "total Tn • 

Physically this quantity represents the number of ions seen by a photon in 
time t n . The parameter Y=l + g 2 /9i is related to tne multiplicities of the 
electronic states g£, for Titanium:Sapphire y=5/3. In this last form of the 
equations we have introduced the possibility of stimulating emission with an 
injection beam I(t). The notation (1-x-Yy) means 

1-x-Yy if l-x-Yy>0 
(l-x-Yy) + = 0 , if l-x-yy< 0 • 

QUALITATIVE ANALYSIS OF THE EQUATIONS 

The purpose of this section is to establish a number of qualitative 
properties of the solutions to the system [10]. We assune the physical 
parameters satisfy 


0 < T 1 <<T l <T ft 

0 < Spo 


Theorem ; If Wp(t), I(t) are positive, continuous and integrable on 
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[0,-hc) and if x(0)>0, y(0) >0, z(0)>0 then 

(i) x(t) >0, y(t)>0, z(t)>0 for all t; 

(ii) x(t), y(t), z(t) are integrable on (0, + -»); 

(iii) The system (10) is asymptotically staole; 

(iv) If x(0)+yy(o) <1, then x(t)+yy(t)<l for all t; and 

(v) If x(0)+yy(o) >1, then there exits T 

such that x(t)+yy(t)<l for all t > T. 

Proof ; We establish the parts of the theroem in the order (iv), (v), (i), 
(ii), (iii). 

(vi) Since x(0)+yy(0)<l, then, by continuity, x(t)+yy(t)<l for t close 
to 0. Suppose x(T)+yy(T)=l and x(t)+yy(t)<l for t<T. Then, at 
t=T, 


(x+yy) (c^+ "^^ 2 ^ )^'*’(c^2"*’"^ < '22 ^ * 

Since <o and C ^ + y C 22 <0, then, if x(T) > 0 and y (T)>0 and not 

both are zero, we would have 

— (x+yy) <0 
dt 

at t=T which would contradict x(t)+yy(t) <1 for t<T. Since x(T)+yy(T)=l, 
not both x(T) and y(T) are zero. 

Suppose y(T)<0. Since y(0)>0, there exists T <T such that y(T^)= 0 
and y(t)>0 for OctcT^. Then, at t=T^, 
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dy 

dt 


= C21 X + axz = (c 91 +az)x. 


21 


Suppose x(T^) < 0 and + az(T^) >0. Then there exists T 2 < T^ 
that x (T^) = 0 and x(t) > 0 for 0 < t < T 2> Then at t = T^ 

^ = c 12 y(T 2 ) + W p (T 2 ) (l-yy(T 2 )) 


such 


dt 


> 0 


since c^ > 0, y(T 2 ) > 0, w p (T 2 ) > 0 and 1 - yy(t 2 ) > 0. 

Hence x increases to the right of T 2> It follows that x cannot be negative 
for 0 < t < T^. So x(T^) > 0. 

Now suppose c^ + az(T^) < 0. Since c^ > 0, a > 0, then z(T^) < 0. 

So there exists T^ < T^ such that z( T^ ) = 0 and z(t) > 0 for 0 < t < T^. 
Then, at t = T^, 


C 31 x + c 32 " y + > 0 

dt 


since I(T 3 ) > 0 and all other quantities are non-negative. Hence z 
increases to the right of T^. So it follows that z(T^) > 0. 


Now suppose x(T^) = 0 so that 


— = 0 at t = T, . Then, at t = T, , 
dt 1 1 


— = W (TO > 0. 
dt p 1 


Additionally, at t = T^, 
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dx . „ dy . ( dx , „ dz^ 
pi — + ^po — + a(z — + x — ) 

dt u at at dt 

n H p ( Tl ) ♦ azUjliyT^X) 


since c 21 >0, W p (T 1 )>0, a>0 and zlT^O. Thus y increases to the right of 
of It now follows that y(T) >0. 

Now suppose x(T)<0. Choose < T such that x(T^) = 0 and x(t)>0 for 
0>t>T^. Then, at t=T^ , 

— = c ?1 >+W (TJd-yy) >0 
dt p 


since y(T 1 )>0, c 21 >0, W ( T ^ ) >0 and l-yy(T 1 )>0. Hence x increases to the 

right of T^. It follows that x (T)>0. 

Since at least one of x(T), y(T) is strictly positive, then at t=T , 

— (x+yy)<0. 
dt 

So x+yy is decreasing at t=T wnich contradicts x+yy<l for t<T . Tnus 
x(t)+yy(t) <1 for all t. 

(v) The proof of (iv) snows tnat if x(0)+yy(G)=l, then x(t)+yy(t)<l 
for all t >0. 

Suppose x(0)+yy(0)>l. Then, for t close to 0, 


d_ 

dt 


(x+yy) =(c 11 +yc 21 )x+(c 12 +yc 22 )y 
<0 


by an argument similar to that given in (i). Thus x+yy is decreasing for t 


close to 0. 

Suppose there exists X >1 such tnat x+yy>x>l for all t. Then 


ii 


3t ( x+y y) = (c 11 +YC 21 )(x+Yy) 

* (c 12 nc 22- 1,(c u nc 21 ))j ' 

< (c^+YC^)^ 

since C 12 +YC 22" Y ^ C 11 +YC 21^ <0 and C 11 + YC 21 <0 ' But then 


x(t)-hry(t)-(x(0)+Yy(0))<(c 11 +yc 21 )At. 


It follows that 


ll£(x(t)+Yy(t) * - * 

which contradicts x +Yy>x >1 for all t. 

Tnus there exists T such tnat x(T) +Yy(T)=l and x(t)+Yy(t)<l for t 
close to T. The analysis in (i) now applies to show that x(t)+Yy(t)<l for 
all t>T. 

(i) The arguments in (iv) show that x,y, are non-negative on an 

interval where x+Yy<l. Hence, if x(0)+Yy(0) <1, then x, y, are 
non-negative for all t. Suppose there exists T such that z(T)<0. 
Choose T^<T such that z(T^)=0. Then, at t=T^, 

c 31 x+c 32 y * I(T l ,>0 

at 

since c 3 j>0, c 32 >0,x>0,y>0 and I (T^ ) >0 . 

Hence z increases to the right of T^ contradicting z(T)<0. Thus, 
z(t)>0 for all t. 

To complete the picture, suppose x(0)+Yy(0)>l. Choose T such tnat 
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x(T)+yY(T)= 1 and x(t)+yy(t)>l for 0<t<T. Then suppose x(T 1 )=0, T^T, and 
x(t)>0 for CKt<T^. Then at t=T^, 


dx 

* ‘ Cl2y ' 

Since c^ 2 >0 and y(T^)>0, x increases to the right of T^. Hence x>0 for all 
t. 

Now suppose y(T 2 )»0, T 2 <T and y(t)>0 for 0 < t <T r At t = T 2 , 
dy 

— = C 21 x + axz 
dt 

Again c^ >0, x>0 and z>0. So y increases to the right of T^. Hence y>, 0 
for all t. 

(ii) Fix T so that x(t)+yy(t)<l for t>T. Then 

— (x+yy) = (c n +YC ?1 )(x+yy) 
dt 11 dL 

+ 

^ c 12 +YC 22" Y ^ c ll +YG 21^ y 

+ 

W (1-x-yy) 

< ( c n + YC 21 )(x + Yy) 

+ 

W (1-x-yy) 

since c^+yc 22 -y(c 22 +yc 2 ^)<0 and y>0. 
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Hence, letting 

E(t) = exp(/ T t (^l p -c 11 -YC 21 )dr) , 

we have 

— (E(t) (x+yy) )< E(t) W . 
dt p 

Thus, 

E(t) (x+yy) - fx(T)+yy(T))</j E(u)W p (u)du. 

It follows that 

0<x+yy<(E(t) _1 [/j E(u) W p (u)du + x(T)+yy(T)]. 

(E(t)-E(u) -1 wJ p (u) du + (E(t) )“ 1 (x(T)+Yy(T)). 
= R(t,T) 

Hence 

0</| (x(t)+yy(t) )dt</j K(t,T)dt. 

Note that, since W > 0, 

P 

exp(-/* W p (r)dr)<l. 

Thus 

0</j (x(t)+yy(t) )dt 


<Jj W p (u)Jy exp((c 1;l +YC 21 )(t-u) )dtdu 


+ 

/j exp((c 11 +yc 21 ) (t-T)) (x(T)+yy(T) )dt 


= / T W p (u) 


exp((c n +yc 21 )(s-u))-l 

c ll +yc 21 


du 


(x(T)ny(T)) 


exp((c 13 +yc 


c 11 +yc 


21 ) (s-T) )-l 
21 


Since + yc^ ^ dnci W p is inte 9 rab1e on [0,+*]» it follows that x+yy is 
integrable on [0,+*]. Consequently both x and y are integrable since each 
is non-negative. 

To see that z is integrable, note that 


z(t) - exp(/* (c 33 +ax) dr)z(T) 

* exp(/* (c 33 +ax)dr) /*exp(-/“(c 33 +ax)dr) (c 31 x+c 32 y+I)du 

Using the fact that x, y, I are non-negative and integrable, an argument 
similar to the one above snows that z is integrable. 


(iii) Let X , Y , Z solve the system (10). Let x=X+ X , 

' ooo o 

y=Y+y Q , z=Z+z q . Assuming x o (0)+yy o (0) <1 and x(0)+yy(0) <1, we 
ODtain tne system 
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f * ( c 11-» yz o- H p>X + < c 12' tH p ) 3C 

-ayx Q Z-ayXZ 

H V 

— = (c21 +a ^o)- + ^22 axZ 0 +a - Z 

dt <- 

= ( c 3i + °ZQ)X+C22Y + (Co3 + aXg)Z+oXZ 
dt 

with X=U=Y_=Z as a solution. 

The linear system 


dX 

dt 


'11 


X + c 



— = c 9 ,Y + c 99 Y 
dt dL ~ 


dZ 

dt 


c 31- + C 32- C 33 Z 


is asymptotical ly stable since all eigenvalues of the coefficient matrix 
have negative real parts. 

Since x q , y Q , z q , W p are non-negative and integrable, it follows from 
Theorem 4.4 of [7 J that X=0=Y=Z is an asymptotically stable solution of tne 
non-linear system. 


NUMERICAL ANALYSIS 

In this section, we discuss the numerical solution to equations (5), 
(7), (8) and (9) subject to the quiescent initial conditions. 

The equations were solved numerically on a DEC VAX 11/750 using a 
Runge-Kutta-Fehlberg (RKF) algorithm [8] and also on a CYBER CY173 using a 
backward difference method) [9]. The time interval over which the solutions 
were computed was from 0 to 200 ns. Material parameters for Ti:Sapphire 
used in the examples are given in Table 1. We were interested in the 
qualitative behavior of the numerical solutions for various punping rates 
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Table 1. Numerical values of the parameters used in these 

calculations. Values come from measurements reported 
in Refs. 5 and 6. 


PARAMETERS 

02=2 

01=3 

SIGMA=3.0E-19 

INDEX OF REFRACTIONS. 76 

NTOTAL=1.0E18 

TAUFL=3 800ns 

TAUl=1.0ns 

TAUC=16.5ns 

TAU2=3000ns 
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and injection strengths as well as the efficiency of tne numerical 
procedures. 

In the examples described here the pimp and injection signals were ooth 
assumed to be Gaussian beams given by 

W (t)= W f(t,t , x ) and 
P P P 

I(t) = I f(t,t 1 ,t 2 j where 
l r (t,t 0 ,x)=2(ln2/ir) 1/2 exp(-ln2((t-t 0 )/.5*T) 2 ) 

with t =30 ns, t =10 ns, t =60 ns, and t =20ns. The values of W and I were 
p p 1 2 

varied. Increased values of W and I correspond to increases in the energies 

delivered d y the punp and injection sources, respectively. 

For relatively low values of W, the RKF algorithm was used to compute 

tne solutions. All of the computations were done in douole precision. In 

-8 

the examples, an error tolerance of 10 was used, but the qualitative 

-5 

behavior of the solutions was preserved for tolerances as large as 10 . 

Figure 2 shows a typical result of the computation. 

The normalized population inversion and photon density are plotted as a 
function of time for a system in which the pumping rate W has the value .01 
ns 1 and I has the value 0, i.e., no injection signal is present. 

Initially, the number of ions in the upper energy levels and the photon 
density in tne medium is zero. At t=0, the pump is turned on and population 
inversion increases, following the leading edge of the pump pulse, until it 
has aosorbed the photons in the pump beam and it levels off. During tnis 
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time the photon density has been increasing steadily (but not visibly on the 
graph). The dynamics up to this point are dominated by the linear part of 
the rate equations. The photon density increases until tne non-linear term 
(tt 4>) becomes comparable to the linear terms and it begins to dominate the 
dynamics. When this occurs there is a rapid increase in the photon density 
accompanied by a corresponding decrease in the population density. This 
latter exchange of photons for excited states oy stimulated emission forms 
the shape of the laser pulse. 

We varied the pumping rate W from .01 ns * to 1.0 ns*'' - by orders of 
magnitude. It is evident from Fig. 3 that the larger the punping rate, the 
earlier stimulated emission will occur and the stronger it will be. The RKF 
algorithm was found to perform very efficiently in the case where W=.01 ns * 
requiring only 822 evaluations of the derivatives. However, as the pump 
rate was increased, the number of derivative evaluations became much larger, 
with 1944 and 6612 evaluations needed when W=.l ns’* and W=1.0 ns"*, 
respectively. The growth of the nunber of derivative evaluations with pump 
rate is shown in Figure 5. 

For large values of W, the RKF algorithm failed to compute the solu- 
tions over the entire time interval in a reasonable number of iterations. W 
was set to 500 ns’* and the RKF procedure required in excess of 67,000 
derivative evaluations to compute the solutions from 0 to 44 ns. Using a 
backward difference method, the solutions, shown in Fig. 4, were found over 
the entire time interval 0 to 200 ns but 17,000 evaluations of the 
derivatives and 3700 evaluations of the Jacobean matrix were necessary. 

When the pump rates are high, the system of differential equations becomes 
stiff, resulting in considerable effort necessary to compute tne solutions. 

In order to examine the effect of the injection beam on the onset of 
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lasing we systematically varied the strength of the injection pulse. Figure 
6 shows the system with W held at .01 ns ^ but varying the value of I from 0 
-1 _! 

to I=1.0e-06 ns and I=1.0e-04 ns . Clearly, the presence of the 
injection beam speeds up the onset of stimulated emission. The RKF 
algorithm worked efficiently in these examples, with each case requiring no 
more than 846 derivative evaluations. 

The theorem of the previous section demonstrates that this system of 
equations is asymptotically stable and the numerical calculations evidence 
this same stability for a wide range of values of tne pimp intensity. The 
non-linearity of the system will support relaxation oscillations at the 
onset of lasing before the stable state is reached. An illustration of these 
oscillations is shown in Fig. 7. The variations of the inversion and photon 
density are shown in separate graphs and also combined together in a pnase 
portrait. An extensive study of the relaxation oscillations of this system 
will appear elsewhere [10]. 


SUMMARY AND CONCLUSION 

We have developed a model of the dynamics of an end-pumped injection 
seeded, four level laser. The model was developed in order to examine the 
early transient behavior of tunable solid state lasers and the effects of 
injection seeding on the timing and shape of the output laser pulse. In 
this paper we report on our analysis of the quantitative behavior of 
solution to the model and on the quality of the nunerical computation of 
these solutions. 

We nave demonstrated tne existence of integraDle solutions to tne 
system of equations forming the model provided that the pump and injection 
source functions are sufficiently regular and that tne initial conditions 
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are physically realizable. Furthermore, we show that these solutions are 
asymptotically stable and establish bounds on them. 

The equations were solved numerically and the sensitivity of the 
numerical calculations to input parameter values was studied. A general 
tendency toward stiffening of the system with increasing pump energy was 
noted. The validity of the model in the anticipated operating regime was 
established. 

ADDENDUM 

The material in this report is presently being prepared for submission 
to the journal. Mathematical Modeling . Co-authors are A. M. Buoncristiani 
and L. F. Roberts. 
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